Measurement of the atmospheric neutrino energy spectrum from 100 GeV to 400 TeV 

with IceCube 



R. Abbasi, 28 Y. Abdou, 22 T. Abu-Zayyad, 33 J. Adams, 16 J. A. Aguilar, 28 M. Ahlers, 32 K. Andeen, 28 J. Auffenberg, 39 
X. Bai, 31 M. Baker, 28 S. W. Barwick, 24 R. Bay, 7 J. L. Bazo Alba, 40 K. Beattie, 8 J. J. Beatty, 18 - 19 S. Bechet, 13 
J. K. Becker, 10 K.-H. Becker, 39 M. L. Benabderrahmane, 40 S. BenZvi, 28 J. Berdermann, 40 P. Berghaus, 28 
D. Berley, 17 E. Bernardini, 40 D. Bertrand, 13 D. Z. Besson, 26 M. Bissok, 1 E. Blaufuss, 17 J. Blumenthal, 1 
D. J. Boersma, 1 C. Bohm, 34 D. Bose, 14 S. Boser, 11 O. Botner, 37 J. Braun, 28 S. Buitink, 8 M. Carson, 22 D. Chirkin, 28 
B. Christy, 17 J. Clem, 31 F. Clevermann, 20 S. Cohen, 25 C. Colnard, 23 D. F. Cowen, 36 - 35 M. V. D'Agostino, 7 
M. Danninger, 34 J. C. Davis, 18 C. De Clercq, 14 L. Demirors, 25 O. Depaepe, 14 F. Descamps, 22 P. Desiati, 28 

G. dc Vries-Uiterweerd, 22 T. DeYoung, 36 J. C. Dfaz-Velez, 28 M. Dierckxsens, 13 J. Dreyer, 10 J. P. Dumm, 28 
M. R. Duvoort, 38 R. Ehrlich, 17 J. Eisch, 28 R. W. Ellsworth, 17 O. Engdegard, 37 S. Euler, 1 P. A. Evenson, 31 

O. Fadiran 4 A. R. Fazely, 6 A. Fedynitch, 10 T. Feusels, 22 K. Filimonov, 7 C. Finley, 34 M. M. Foerster, 36 B. D. Fox, 36 
A. Franckowiak, 11 R. Franke, 40 T. K. Gaisser, 31 J. Gallagher, 27 M. Geisler, 1 L. Gerhardt, 8 ' 7 L. Gladstone, 28 
T. Gliisenkamp, 1 A. Goldschmidt, 8 J. A. Goodman, 17 D. Grant, 21 T. Griesel, 29 A. GroB, 16 ' 23 S. Grullon, 28 
M. Gurtner, 39 C. Ha, 36 A. Hallgren, 37 F. Halzen, 28 K. Han, 16 K. Hanson, 13 ' 28 K. Helbing, 39 P. Herquet, 30 
S. Hickford, 16 G. C. Hill, 28 K. D. Hoffman, 17 A. Homeier, 11 K. Hoshina, 28 D. Hubert, 14 W. Huelsnitz, 17 'E 
J.-P. Hiilfi, 1 P. O. Hulth, 34 K. Hultqvist, 34 S. Hussain, 31 A. Ishihara, 15 J. Jacobsen, 28 G. S. Japaridze, 4 
H. Johansson, 34 J. M. Joseph, 8 K.-H. Kampert, 39 A. Kappes, 9 T. Karg, 39 A. Karle, 28 J. L. Kelley, 28 N. Kemming, 9 

P. Kenny, 26 J. Kiryluk, 8 ' 7 F. Kislat, 40 S. R. Klein, 8 ' 7 J.-H. Kohne, 20 G. Kohnen, 30 H. Kolanoski, 9 L. Kopke, 29 
D. J. Koskinen, 36 M. Kowalski, 11 T. Kowarik, 29 M. Krasberg, 28 T. Krings, 1 G. Kroll, 29 K. Kuehn, 18 T. Kuwabara, 31 
M. Labare, 14 S. Lafebre, 36 K. Laihem, 1 H. Landsman, 28 M. J. Larson, 36 R. Lauer, 40 R. Lehmann, 9 J. Liinemann, 29 
J. Madsen, 33 P. Majumdar 40 A. Marotta, 13 R. Maruyama, 28 K. Mase, 15 H. S. Matis, 8 M. Matusik, 39 K. Meagher, 17 
M. Merck, 28 P. Meszaros, 35 - 36 T. Meures, 1 E. Middell, 40 N. Milke, 20 J. Miller, 37 T. Montaruli, 28 '[l] R. Morse, 28 
S. M. Movit, 35 R. Nahnhauer, 40 J. W. Nam, 24 U. Naumann, 39 P. Niefien, 31 D. R. Nygren, 8 S. Odrowski, 23 
A. Olivas, 17 M. Olivo, 37 ' 10 A. O'Murchadha, 28 M. Ono, 15 S. Panknin, 11 L. Paul, 1 C. Perez de los Heros, 37 
J. Pctrovic, 13 A. Piegsa, 29 D. Pieloth, 20 R. Porrata, 7 J. Posselt, 39 P. B. Price, 7 M. Prikockis, 36 G. T. Przybylski, 8 
K. Rawlins, 3 P. Redl, 17 E. Resconi, 23 W. Rhode, 20 M. Ribordy, 25 A. Rizzo, 14 J. P. Rodrigues, 28 P. Roth, 17 
F. Rothmaier, 29 C. Rott, 18 T. Ruhe, 20 D. Rutledge, 36 B. Ruzybayev, 31 D. Ryckbosch, 22 H.-G. Sander, 29 
M. Santander, 28 S. Sarkar, 32 K. Schatto, 29 S. Schlenstedt, 40 T. Schmidt, 17 A. Schukraft, 1 A. Schultes, 39 O. Schulz, 23 
M. Schunck, 1 D. Seckel, 31 B. Semburg, 39 S. H. Seo, 34 Y. Sestayo, 23 S. Seunarine, 12 A. Silvestri, 24 K. Singh, 14 
A. Slipak, 36 G. M. Spiczak, 33 C. Spiering 40 M. Stamatikos, 18 <|l] T. Stanev, 31 G. Stephens, 36 T. Stezelberger, 8 
R. G. Stokstad, 8 S. Stoyanov, 31 E. A. Strahler, 14 T. Straszheim, 17 G. W. Sullivan, 17 Q. Swillens, 13 H. Taavola, 37 
I. Taboada, 5 A. Tamburro, 33 O. Tarasova 40 A. Tepe, 5 S. Ter-Antonyan, 6 S. Tilav, 31 P. A. Toale, 36 S. Toscano, 28 
D. Tosi, 40 D. Turcan, 17 N. van Eijndhoven, 14 J. Vandenbroucke, 7 A. Van Overloop, 22 J. van Santcn, 28 M. Voge, 23 
B. Voigt, 40 C. Walck, 34 T. Waldenmaier, 9 M. Wallraff, 1 M. Walter, 40 Ch. Weaver, 28 C. Wendt, 28 S. Westerhoff, 28 
N. Whitehorn, 28 K. Wiebe, 29 C. H. Wiebusch, 1 G. Wikstrom, 34 D. R. Williams, 2 R. Wischnewski, 40 

H. Wissing, 17 M. Wolf, 23 K. Woschnagg, 7 C. Xu, 31 X. W. Xu, 6 G. Yodh, 24 S. Yoshida, 15 and P. Zarzhitsky 2 

(IceCube Collaboration) 

III. Physikalisches InsUtut, RWTH Aachen University, D-52056 Aachen, Germany 
2 Department of Physics and Astronomy, University of Alabama, Tuscaloosa, Alabama 35487, USA 
3 Department of Physics and Astronomy, University of Alaska Anchorage, 
3211 Providence Dr., Anchorage, Alaska 99508, USA 
4 CTSPS, Clark-Atlanta University, Atlanta, Georgia 30314, USA 
5 School of Physics and Center for Relativistic Astrophysics, 
Georgia Institute of Technology, Atlanta, Georgia 30332, USA 
6 Department of Physics, Southern University, Baton Rouge, Louisiana 70813, USA 
7 Department of Physics, University of California, Berkeley, California 94720, USA 
8 Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA 
9 Institut fur Physik, Humboldt- Universitat zu Berlin, D- 12489 Berlin, Germany 
10 Fakultdt fur Physik & Astronomie, Ruhr- Universitat Bochum, D-44 780 Bochum, Germany 
11 Physikalisches Institut, Universitat Bonn, Nussallee 12, D-53115 Bonn, Germany 
12 Department of Physics, University of the West Indies, 
Cave Hill Campus, Bridgetown BB11000, Barbados 
13 Universite Libre de Bruxelles, Science Faculty CP230, B-1050 Brussels, Belgium 
1J> Vrije Universiteit Brussel, Dienst ELEM, B-1050 Brussels, Belgium 



2 



15 Department of Physics, Chiba University, Chiba 263-8522, Japan 
16 Department of Physics and Astronomy, University of Canterbury, Private Bag 4800, Christchurch, New Zealand 
"Department of Physics, University of Maryland, College Park, Maryland 20742, USA 
18 Department of Physics and Center for Cosmology and Astro-Particle Physics, 

Ohio State University, Columbus, Ohio 43210, USA 
19 Department of Astronomy, Ohio State University, Columbus, Ohio 43210, USA 
20 Department of Physics, TU Dortmund University, D-44^21 Dortmund, Germany 
21 Department of Physics, University of Alberta, Edmonton, Alberta, Canada T6G 2G7 
22 Department of Subatomic and Radiation Physics, University of Gent, B-9000 Gent, Belgium 
28 Max- Planck- Institut fur Kernphysik, D-691 77 Heidelberg, Germany 
2/> Department of Physics and Astronomy, University of California, Irvine, California 92697, USA 
25 Laboratory for High Energy Physics, Ecole Polytechnique Federale, CH-1015 Lausanne, Switzerland 
26 Department of Physics and Astronomy, University of Kansas, Lawrence, Kansas 66045, USA 
27 Department of Astronomy, University of Wisconsin, Madison, Wisconsin 53706, USA 

28 Department of Physics, University of Wisconsin, Madison, Wisconsin 53706, USA 
29 Institute of Physics, University of Mainz, Staudmger Weg 7, D-55099 Mainz, Germany 
80 Universite de Mons, 7000 Mons, Belgium 
Bartol Research Institute and Department of Physics and Astronomy, 
University of Delaware, Newark, Delaware 19716, USA 
32 Department of Physics, University of Oxford, 1 Keble Road, Oxford 0X1 3NP, UK 
33 Department of Physics, University of Wisconsin, River Falls, Wisconsin 54022, USA 
SJt Oskar Klein Centre and Department of Physics, 
Stockholm University, SE- 10691 Stockholm, Sweden 
Department of Astronomy and Astrophysics, Pennsylvania State University, University Park, Pennsylvania 16802, USA 
86 Department of Physics, Pennsylvania State University, University Park, Pennsylvania 16802, USA 
"Department of Physics and Astronomy, Uppsala University, Box 516, S-75120 Uppsala, Sweden 
88 Department of Physics and Astronomy, Utrecht University /SRON, NL-3584 CC Utrecht, The Netherlands 
89 Department of Physics, University of Wuppertal, D-42119 Wuppertal, Germany 
4 "DESY, D-15735 Zeuthen, Germany 
(Dated: December 20, 2010) 

A measurement of the atmospheric muon neutrino energy spectrum from 100 GeV to 400 TeV 
was performed using a data sample of about 18,000 up-going atmospheric muon neutrino events in 
IceCube. Boosted decision trees were used for event selection to reject mis-reconstructed atmospheric 
muons and obtain a sample of up-going muon neutrino events. Background contamination in the 
final event sample is less than 1%. This is the first measurement of atmospheric neutrinos up to 400 
TeV, and is fundamental to understanding the impact of this neutrino background on astrophysical 
neutrino observations with IceCube. The measured spectrum is consistent with predictions for the 
atmospheric f M + flux. 



PACS numbers: 95.55.Vj,95.85.Ry,14.60.Lm,29.40.Ka 
I. INTRODUCTION 

The IceCube neutrino telescope [l|, currently under 
construction in the glacial ice at the South Pole, is capa- 
ble of detecting high energy neutrinos of all three flavors. 
In particular, charged current (CC) interactions between 
or D^, and nucleons in the ice, produce muons. Ice- 
Cube detects the Cherenkov radiation produced as these 
muons propagate and undergo radiative losses. By recon- 
structing the muon's track and energy loss, the direction 
and energy of the incident neutrino can be inferred. 

Atmospheric neutrinos are produced in the decay 
chains of particles created by the interaction of cosmic 
rays with the Earth's atmosphere [2|t4J- IceCube has 
an unprecedented high statistics, high energy reach for 
these atmospheric neutrinos. Hence, IceCube can be 
used to test predictions for the flux of atmospheric neutri- 
nos at high energies, including the uncertain contribution 
from charm production above about 100 TeV. The atmo- 



spheric neutrino flux can also be used to verify that the 
IceCube detector is performing as expected @- Under- 
standing the energy and zenith dependence of the atmo- 
spheric neutrino flux in IceCube is important since this is 
an irreducible background for searches for a diffuse flux, 
or for point sources, of astrophysical neutrinos. 

This analysis used data taken from April 2008, to May 
2009, while IceCube was operating in a 40-string con- 
figuration. The signal events were up-going atmospheric 

and interactions. The background was down-going 
atmospheric muons that were mis-reconstructed as up- 
going. An as- yet unmeasured but anticipated diffuse flux 
of astrophysical neutrinos was ignored. Predictions for 
this flux are negligible compared to predictions for the 
atmospheric neutrino flux, over most of the energy range 
for this analysis, and it can readily be accommodated 
within the reported uncertainties. 

Boosted decision trees (BDT) were used to obtain an 
event sample with negligible background contamination. 



An unfolding of the atmospheric neutrino energy spec- 
trum, over the neutrino energy range 100 GeV to 400 
TeV, was performed. Systematic uncertainties in the un- 
folded spectrum were estimated and highlight the efforts 
that are underway to reduce systematic uncertainties in 
neutrino measurements with IceCube. 

We will briefly review the production and distribution 
of atmospheric neutrinos in Section [II] In Section IIII1 
we will discuss the IceCube detector, and the detection 
of muon neutrino events in IceCube. Event reconstruc- 
tion and event selection specific to this analysis will be 
discussed in Section HVl The unfolding analysis and sys- 
tematic uncertainties will be discussed in Section [V[ Fi- 
nally, in Section I VII we will discuss the implications of 
the unfolded result. 



II. ATMOSPHERIC NEUTRINOS 

Cosmic rays are high energy particles, mostly protons 
and helium nuclei, but also heavier ionized nuclei, that 
are believed to be accelerated in various astrophysical 
phenomena @, 0] ■ Possible cosmic ray production sites 
include active galactic nuclei, gamma ray bursts, and 
supernova explosions. Detecting astrophysical neutri- 
nos, produced in conjunction with cosmic rays at point 
sources such as these, is one of the primary goals of Ice- 
Cube. The energy spectrum of cosmic rays is rather 
steep, dN/dE oc E~ 2 - 7 , and steepens to dN/dE oc E~ 3 
above the "knee", or about 10 6 GeV Q. A possible sec- 
ond knee is a steepening to about E~ 3 - 2 above 5 x 10 8 
GeV [8]. A further kink in the spectrum has been ob- 
served at ~ 3 x 10 9 GeV, where the spectrum flattens to 
dN/dE oc E~ 2 ' 7 again. The event sample for this anal- 
ysis is primarily the result of interactions of cosmic rays 
with energies below the first knee. 

Hadronic interactions between cosmic rays and parti- 
cles in the Earth's atmosphere produce large numbers 
of mesons, primarily pions and kaons. Hundreds or even 
thousands of these mesons can be produced in the shower 
that follows the interaction of a single high energy cos- 
mic ray. Neutrinos are produced in the leptonic or semi- 
leptonic decays of charged pions or kaons, as well as in 
the subsequent decay of the muons. Neutrinos from muon 
decay are important up to a few GeV. Pions and kaons 
that decay in-flight are the primary source of atmospheric 
muon neutrinos from a few GeV up to about 100 TeV. 
With rest-frame lifetimes on the order of 10 -8 s, these 
mesons often lose some of their energy in collisions prior 
to decaying, leading to lower energy neutrinos among the 
decay products. Hence, the spectral slope of this "con- 
ventional" atmospheric neutrino flux 0, H| asymptoti- 
cally becomes one power steeper than that of the pri- 
mary cosmic ray spectrum. Theoretical uncertainties in 
predictions for the conventional flux are dominated by 
uncertainties in the normalization and spectral distribu- 
tion of the cosmic ray flux. Additional uncertainties in- 
clude the ratio of pions to kaons produced by cosmic ray 
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FIG. 1: The predicted flux of atmospheric muon neutrinos. 
The solid line is the conventional + flux [3j], averaged 
over the zenith range 90° to 180°. The long-dashed line is 
the prompt + flux [Sj, also averaged over the zenith 
range 90° to 180°. The dot-dashed curve is the sum of the 
conventional and prompt models. The flux predictions from 
[j| were extended to higher energies as discussed in Sect. lIII Cl 
For reasons discussed in Sect. IV Bl the zenith region from 90° 
to 97° was not used in the analysis. The zenith-averaged 
conventional flux, for the range 97° to 180°, is included in 
the figure as the small-dashed line. The prediction for the 
zenith-averaged prompt flux is not affected by this change in 
angular region. 



interactions, which affects the zenith angle distribution, 
particularly near the horizon. 

At sufficiently high energies, another production mech- 
anism is possible. The "prompt" atmospheric neutrino 
flux [9l4llj is made up of neutrinos produced in the semi- 
lcptonic decays of charmed mesons and baryons. These 
particles decay almost immediately (rest-frame lifetimes 
on the order of 10~ 12 s), before losing energy in colli- 
sions. Hence, the spectrum for the prompt flux more 
closely follows the cosmic ray spectrum and is about one 
power harder than the conventional flux at high energy. 
The prompt flux has not yet been measured, but is ex- 
pected to be important above about 100 TeV @,[IU. Just 
like the conventional flux, predictions for the prompt flux 
are impacted by uncertainties in the normalization and 
spectral distribution of the cosmic ray flux. Additional 
sources of uncertainty for t he p rompt flux include charm 
production cross sections [131 ] and fragmentation func- 
tions, which have not been measured at these energies 
in accelerator experiments. Figure [T] shows the predicted 
flux of conventional and prompt atmospheric muon neu- 
trinos [3J, |9| . 

Although high energy cosmic rays arrive almost 
isotropically, with deviations less than 0.1% [lij], the 
zenith angle dependence of high energy atmospheric neu- 
trino production is complicated by the direction of the 
shower through the atmosphere. The energy spectrum 
of nearly horizontal conventional atmospheric neutrinos 
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is flatter than that of almost vertical neutrinos because 
pions and kaons in inclined showers spend more time in 
the tenuous atmosphere where they are more likely to 
decay before losing energy in collisions. Additionally, at- 
tenuation of the neutrino flux by the Earth is a function 
of energy and zenith angle. Above about 10 TeV, attenu- 
ation of the neutrino flux in the Earth is important, and 
affects the zenith and energy dependence of the flux at 
the detector. 



III. NEUTRINO DETECTION WITH ICECUBE 

A. The IceCube Detector 

IceCube p], is able to detect neutrinos over a wide 
energy range, from about 100 GeV to more than 10 9 GeV. 
The design is a balance between energy resolution, angu- 
lar resolution, energy range, and cost, and was driven by 
the goal of detecting astrophysical neutrino point sources, 
which are believed to be correlated with cosmic ray pro- 
duction sites. A large detector is required as a result of 
the extremely small cross-sections for neutrino interac- 
tions, as well as the low fluxes expected for astrophysical 
neutrinos. 

When completed in 2011, IceCube will comprise 86 
strings, with 5160 photomultiplier tubes (PMT). Each 
string includes sixty digital optical modules (DOM). A 
DOM is a single PMT and associated electronics in a 
glass pressure sphere. The instrumented part of the array 
extends from 1450 m to 2450 m below the surface of the 
ice. Horizontally, 78 of the strings are 125 m apart and 
spread out in a triangular grid over a square kilometer, 
so that the entire instrumented volume will be 1 km 3 of 
ice. Vertical DOM spacing is a uniform 17 m for these 78 
strings. A subset of the detector, known as "DeepCore" , 
consists of eight specialized and closely spaced strings of 
sensors located around the center IceCube string. 

Figure [5] shows the IceCube observatory and its com- 
ponent arrays. This analysis used data from 359 days 
of livetime while operating in a 40-string configuration, 
from April 2008, to May 2009. Figure [3] shows an over- 
head view of the layout of the 40-string configuration, 
which was roughly twice as long in one horizontal direc- 
tion as in the other. 

At the heart of each DOM is a 10 inch (25 cm) Hama- 
matsu PMT [16| (see Fig.[4j. A single Cherenkov photon 
arriving at a DOM and producing a photoelectron is de- 
fined as a hit. DOM main board electronics [ljj apply a 
threshold trigger to the PMT analog output. This thresh- 
old is equivalent to 0.25 of the signal generated by a pho- 
toelectron, after amplification by the PMT. When this 
threshold is exceeded, local coincidence checks between 
this DOM and nearest neighbor or next-to-nearest neigh- 
bor DOMs on a string are performed to reduce false trig- 
gers that result from dark noise. If a nearest or next-to- 
nearest neighbor DOM also has a detection above thresh- 
old within a ±1000 ns window, the PMT total charge 



waveforms are digitized, time stamped, and sent to the 
surface. The digitized waveform from a DOM can con- 
tain several pulses, and each pulse can be the result of 
multiple photoclectrons. The simple majority trigger for 
building an event is eight hit DOMs within a 5000 ns 
trigger window. 

The data rate from the data acquisition system at 
the South Pole far exceeds the amount of data that can 
be transmitted via satellite. Hence, a significant reduc- 
tion in the trigger-level data must be accomplished with 
software-based filtering at the South Pole. A cluster of 
processors performs a variety of fast reconstructions on 
the data, and applies multiple software-based filters to 
the results. These filters either reject events that are 
uninteresting background events, or extract particular 
classes of events. Events are sent to a buffer if they pass 
one or more of the filters. The transfer of data from this 
buffer over a communications satellite is handled by the 
South Pole Archival and Data Exchange (SPADE) sys- 
tem. 

The deep glacial ice at the South Pole is optically trans- 
parent, making it an ideal medium for a large volume 
Cherenkov detector. The ice sheet is just over 2800 m 
thick and was created over a period of roughly 165,000 
years [18]. It serves multiple roles: a stable platform for 
the DOMs, the target medium for neutrino interactions, 
the propagation and detection medium for Cherenkov 
photons produced by charged particles, and an overbur- 
den for attenuation of down-going atmospheric muons. 
Upward moving particles will have had to result from 
particles that penetrated the Earth and can readily be 
identified as resulting from neutrino interactions. 

Optical properties of the ice are discussed in [l9l - [2l"j |. 
Scattering and absorption of photons in the ice is caused 
by bubbles, dust particles, and crystal defects. Below 
about 1400 m, the ice is essentially free of bubbles, and 
scattering is dominated by dust. Micron-sized dust grains 
were carried as wind-borne aerosols during the periods 
of ice formation, and deposited in the ice. Variations 
with depth are due to the periodic build up of dust that 
resulted from the prevailing atmospheric conditions when 
the layers of ice were being formed. 

The depth and wavelength dependence of scatter- 
ing and absorption as measured in the ice around the 
AMANDA detector is discussed in Ref. [nj. Now sur- 
rounded by IceCube and no longer operating, AMANDA 
was the predecessor and prototype for IceCube. Ice prop- 
erties were extrapolated to lower depths using ice core 
measurements taken at Vostok Station and Dome Fuji in 
Antarctica, then scaled to the location of IceCube using 
an age vs. depth relationship Studies are on-going 
that use LEDs on flasher boards within each DOM to 
directly measure ice properties in the deepest ice instru- 
mented by IceCube. 
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FIG. 2: IceCube Neutrino Observatory and its component arrays. 



B. Muon Neutrino Detection 

Muon neutrinos undergoing CC interactions in the ice 
produce muons. The muons on average carry about 75% 
of the initial neutrino energy [22| . Simulation studies in- 
dicate that muon angular resolution is typically between 
0.5° and 1°, depending on the angle of incidence and the 
muon energy. The energy loss per meter, for a muon 
propagating through the ice, is related to its energy [gj]: 

-{^j=a{E)+P{E)E, (1) 

where E is the muon energy, a « 0.24 GeV/m is the 
ionization energy loss per unit propagation length, and 
f3 ?s 3.3 x 1CP 4 m _1 is the radiative energy loss through 
bremsstrahlung, pair production, and photonuclear scat- 
tering, (a and f3 are both weak functions of energy). For 
muon energies less than about a TeV, energy loss is dom- 
inated by ionization, and the light produced is nearly in- 
dependent of energy. However, for higher energy muons, 
there are many stochastic interactions along the muon's 
path and there is a linear relationship between the en- 



ergy loss per meter and the muon energy. Most of the 
Cherenkov light emitted along the muon's path comes 
from the secondary particles produced in radiative losses. 
An estimation of dE/dX, based on the amount of de- 
tected light, the event geometry, and the ice properties, 
was used in the energy spectrum unfolding discussed in 
Sect. El The energy of individual events was not esti- 
mated. Rather, the distribution of neutrino energies was 
directly inferred from the distribution of reconstructed 
muon dE/dX values. 



The detection rate for high energy (y^) is aided by 
the fact that the CC interaction cross section, as well as 
the range of the resultant muon, are proportional to the 
neutrino energy. High energy muons have a significant 
path-length and can reach the detector even if produced 
outside of the detector, hence increasing the effective vol- 
ume. Muons in earth or ice can have a track length from 
several tens of meters, up to several kilometers, depend- 
ing on the muon energy and the detection threshold. The 
average track length, before the muon energy falls below 
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FIG. 3: Overhead view of IceCube 40 string configuration. 
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a detection threshold , is given by: 



4 in 



a + QEf 



(2) 



where En is the initial muon energy. 



C. Simulation 

Simulation of atmospheric muons and neutrinos was 
used for determining event selection and background re- 
jection cuts. Simulation was also used for the response 
matrix (discussed in Section|V| and the predicted dE/dX 
distribution for the unfolding analysis. Several special- 
ized simulated data sets were used for systematics studies 
and toy Monte Carlo (MC) studies. 



Muons from air showers were simulated with COR- 
SIKA 123] . The primary cosmic ray energy spectrum 
known as the Horandel poly-gonato model [8] was used. 
In this model, the spectrum of each component is a com- 
bination of two power laws, with the turnover between 
the two power laws being a function of the nuclear charge 
Z of the primary cosmic ray. CORSIKA propagates cos- 
mic ray primaries (up to Fe) to their point of interac- 
tion with a nucleus in the atmosphere. Hadronic inter- 
actions in the atmosphere were modelled using the in- 
teraction model SIBYLL [24]]. Secondary particles were 
then tracked until they interacted or decayed. Coincident 
muons in the detector, originating from separate cosmic 
ray events, were accounted for by combining simulated 
events and re-weighting them to account for the proba- 
bility of coincident events occurring. 

Muon propagation and energy loss within and around 
the detector was simulated with the program MMC 
(Muon Monte Carlo) [Hj]. MMC accounts for ioniza- 
tion, bremsstrahlung, photo-nuclear interactions, and 
pair production. In addition to muon tracks and ener- 
gies, secondary particles from the stochastic energy losses 
are included in the output of MMC. The production and 
propagation of Cherenkov light from the muons and sec- 
ondary particles was simulated using the program Pho- 
tonics [26| , which accounts for the depth-dependent scat- 
tering and absorption properties of the ice. Direct track- 
ing of Cherenkov photons through the layered glacial ice 
was too computationally intensive for simulation produc- 
tion. Photonics was run beforehand to create lookup ta- 
bles which were then used during the detector simulation. 
The tables included light yield and photon propagation 
time distributions at a given location in the ice from a 
given source type and location. Simulation of the de- 
tector response to electromagnetic and hadronic showers 
(so-called cascade events) also used pre-tabulated light 
yield tables and photon propagation time information 
generated by Photonics. An energy-dependent scaling 
factor was applied for hadronic cascades, to account for 
the fact that hadronic cascades produce less Cherenkov 
light than their electromagnetic counterparts 27}. 

Neutrino propagation from point of origin in the at- 
mosphere to interaction in or near the detector was sim- 
ulated with ANIS ANIS generates neutrinos of 
any flavor according to a specified flux, propagates them 
through the Earth, and in a final step simulates neutrino 
interactions within a specified volume. All simulated neu- 
trinos were forced to interact, but their probability of in- 
teracting was included in the event weight assigned by 
ANIS. ANIS accounts for CC and neutral current (NC) 
neutrino-nucleon interactions, as well as neutrino regen- 
eration following NC interactions. Also accounted for 
is the offset between neutrino propagation direction and 
the direction of the outgoing muon following a CC inter- 
action. Cross sections for and CC and NC inter- 
actions were based on the CTEQ5 parton distributions 
[29} . The density profile in the Earth was modeled using 
the Preliminary Reference Earth Model [30|. 
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FIG. 5: Effective livetime of the simulation used to train the 
unfolding algorithm, as a function of neutrino energy. 



Simulated neutrino events were generated with an E~ 2 
spectral index, then weighted according to their contri- 
bution to the atmospheric neutrino flux. The flux predic- 
tions of Honda et al. [3j were used for conventional atmo- 
spheric muon neutrinos, and those of Enberg et al. [9( for 
prompt atmospheric muon neutrinos. The predictions for 
muon neutrinos from pions and kaons were extended to 
higher energies by fitting a physics-motivated analytical 
equation based on energy and zenith angle (Ref. [12[ and 
chapter 7 of Ref. 0) in an overlapping region with the 
detailed calculations of Honda et al. [3] . 

Since simulated events were generated with a harder 
spectrum than atmospheric neutrinos, the effective live- 
time for high energy events was boosted. Additionally, 
since all events were forced to interact in or near the de- 
tector, the effective livetime for low energy events was 
boosted. The effective livetime of the neutrino simula- 
tion used to train the unfolding algorithm is shown in 
Fig. El 



IV. EVENT RECONSTRUCTION AND 
BACKGROUND REJECTION 

A variety of algorithms are used in IceCube for event 
reconstruction, classification, and background rejection, 
depending on the energy range, the anticipated signal 
and backgrounds for a particular analysis, as well as 
the neutrino flavor. The background for this analy- 
sis was down-going atmospheric muons that were mis- 
reconstructed as up-going. Despite the depth of IceCube, 
the ratio of down-going atmospheric muons to muons pro- 
duced in or near the detector by neutrino interactions is 
roughly one million to one [l[ • Below is a brief summary 
of the muon track reconstruction algorithms and event 
selection methods that were used in this analysis. 



A. Event Reconstruction 

The LineFit reconstruction is a fast, first-guess algo- 
rithm based on the assumption that the Cherenkov pho- 
tons from a muon propagate on a wavefront perpendic- 
ular to the track. This assumption leads to a fitting al- 
gorithm that is extremely fast, and often estimates the 
muon track direction within ten degrees. LineFit, and 
likelihood-based reconstructions (discussed next) seeded 
with the LineFit track, were used as part of the software 
filtering at the South Pole. Additionally, the wavefront 
velocity estimated by LineFit is correlated with how well 
the track hypothesis fits the distribution of recorded light 
and was used as an event selection cut prior to one of the 
two BDTs. 

Maximum likelihood reconstruction algorithms ac- 
count for the geometric dependence of photon arrival 
times, as well as the stochastic variability in arrival times 
due to scattering in the ice. The likelihood function to 
be maximized is the function [31| 

£ = Ylp( a ,tut,i), (3) 

i 

where a is the set of parameters characterizing the hy- 
pothesized track, i.e. three coordinates for the vertex lo- 
cation, two angles for the direction, and possibly energy, 
and p (a, t^it) is t ne probability distribution function [32[ 
for photon hit times, given the track hypothesis. The 
product is over all photon hits in the event. In practice, 
the maximum of the likelihood function is found by min- 
imizing the negative of the log of the likelihood, so the 
product becomes a sum. To further simplify implemen- 
tation, a transformation is made and time residual, t Ies , 
is used in place of hit time, thit, where 

^rcs — ^hit ^gco- (4) 

The geometric travel time, f geo , is based on a straight 
photon path with no scattering. 

Single Photoelectron (SPE) fits are likelihood recon- 
structions that use only the arrival time of the first pho- 
toelectron in all hit DOMs. Typically, 16 or 32 itera- 
tions of the SPE fit are performed, with the seed track 
randomly altered for each iteration. This helps ensure 
that a local minimum is not chosen as the final track. 
The Multiple Photoelectron (MPE) fit is similar to the 
SPE fit, however, it uses the total number of observed 
photons to describe the arrival time of the first photon. 
When many photons arrive at the same DOM, the first 
photon is scattered less than an average photon. Since 
more information is used, the directional accuracy of the 
fit is often improved slightly, as compared to the SPE 
fit. Moreover, using track quality parameters based on 
the MPE fit rather than on the SPE fit provided better 
event discrimination and improved the signal efficiency 
of the BDTs by about 10%. 

In addition to track location and direction, the like- 
lihood reconstructions return several variables that are 
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FIG. 6: Distribution of the PLogL variable (from the MPE 
fit), for neutrino simulation, muon background simulation, 
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amount of data requiring higher level processing. PLogL was 
then recalculated based on the result of the MPE fit. 




" data 

■ signal (simulation) 

■ bkgnd (simulation) 



3 
ra 

a 



20 
.MPE 



&f \ (degrees) 



FIG. 7: Distribution of the difference between LineFit zenith 
angle and MPE fit zenith angle, for neutrino simulation, muon 
background simulation, and for data. 



used to estimate fit quality. These variables include 
the log-likelihood (LogL) and the reduced log-likelihood 
(RLogL). RLogL = LogL/n<j f, where n<j f is the number 
of degrees of freedom in the minimization, i.e. the num- 
ber of hit DOMs minus the number of parameters to be 
fit. RLogL is then (ideally) independent of the number 
of hit DOMs. A similar scaled parameter called PLogL, 
equal to LogL/(number of hit DOMs - 2.5), has also been 
found to provide additional discriminatory power. Fig- 
ure [6] shows the distribution of the PLogL variable from 
the MPE fit, that was used by the BDTs. 

Figure [7] shows the distribution of the difference be- 
tween LineFit zenith angle and MPE fit zenith angle, for 
signal simulation and for data. This angular difference 
was used as an input to the BDTs. 

PhotoRec is a reconstruction algorithm that accounts 
for spatially variable ice properties [33 |. It does this by 
incorporating light propagation tables created by Pho- 
tonics [26| . The output of PhotoRec used in this analysis 
was the estimation of dE/dX, the average muon energy 



loss per unit propagation length (GeV/m) that would 
produce the detected amount of light. The reconstructed 
dE/dX is proportional to the number of photons de- 
tected and hence to the number of photons emitted along 
the muon's track. To correctly scale the proportionality, 
changes in the photon intensity due to the distance be- 
tween the track and the hit DOMs, and the amount of 
scattering and absorption between light generation and 
detection points, are accounted for. The reconstruction 
algorithm incorporates the detailed ice model, but as- 
sumes that stochastic energy losses are uniform along the 
track. As mentioned in Section llll B[ the amount of light 
emitted along the track of a high energy muon (greater 
than about a TeV) is linearly correlated with the muon's 
energy and it is possible to estimate the muon's energy 
(near the center of the detected track), using dE/dX, 
with an accuracy of 0.3 on a log scale. However, for 
low energy muons, the amount of light emitted along 
the track is nearly independent of energy, as discussed 
in Section IIII B[ Additionally, the PhotoRec algorithm 
does not account for the fraction of detected photons 
that may be from the hadronic shower at the interac- 
tion point (if that occurs inside the detector), nor does 
it account for the length of the muon track inside the de- 
tector. Hence, the correlation between dE/dX and muon 
energy degrades below a TeV. 

The paraboloid algorithm [34j analyzes the value of the 
likelihood function around a seed track. After transform- 
ing the coordinate space to one centered on the direction 
of the seed track, it fits a constant likelihood ellipse to 
the likelihood space around the direction of the track. 
The important result is the paraboloid sigma, calculated 
from the major and minor axes of the constant likeli- 
hood ellipse. Paraboloid sigma provides an estimate of 
the pointing error of the track. 

In a Bayesian reconstruction, the standard likelihood 
function is multiplied by a bias function which depends 
only on the event hypothesis and not on the actual 
event data. The bias is used as a way to include prior 
knowledge of the characteristics of the data, that mis- 
reconstructed down-going tracks dominate the signal by 
about three orders of magnitude at this stage. The 
Bayesian likelihood ratio is the useful result from this re- 
construction, LogL Baycs — LogL gPE32 , where SPE32 refers 
to a 32-iteration SPE fit. Figure |5] shows the distribution 
of the Bayesian likelihood ratio, for simulation and for 
data. 

In a reconstruction algorithm known as the umbrella 
fit, the minimizer is constrained to track directions with 
space angles more than 90° from a seed track. The like- 
lihood ratio LogL Umbl . — LogL SPE3 2 is used as an event 
selection parameter. Good tracks have a higher SPE like- 
lihood than a fit constrained to have a directional com- 
ponent in the opposite direction. This reconstruction 
provides discriminating power for certain events that are 
stuck in a local minimum in the likelihood space, such as 
down-going or near horizontal events that reconstruct as 
directly up-going. 
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FIG. 8: Distribution of the Bayesian likelihood ratio, for neu- 
trino simulation, muon background simulation, and for data. 

Split track reconstructions begin by creating four sub- 
events from the initial event. Two sub-events are created 
by separating all hit DOMs into the group hit before 
the average time, and the group hit after the average 
time. Two additional sub-events are based on geome- 
try. All hit DOMs are projected perpendicularly along 
the track. Then, the DOMs are split into two groups 
based on whether they fall before or after the location 
of the center of gravity of the pulses. LineFit and SPE 
reconstructions are performed on each of these four sub- 
sets. These fits provide discrimination for poorly recon- 
structed tracks, as well as for tracks that reconstruct as 
up-going due to the superposition of hits from two sep- 
arate down-going muons. A loose cut on zenith angles 
from the split track reconstructions was used as an event 
selection cut prior to one of the two BDTs. Addition- 
ally, the zenith angles were used as input variables for 
the BDTs. Figure [5] compares the zenith angles from the 
SPE fits for the two sub-events found by the geometric 
split. 

In addition to zenith angles and likelihood ratios, sev- 
eral other measured or reconstructed variables were used 
for event discrimination. For example, the likelihood that 
a track is properly reconstructed is correlated with the 
number of hit strings (NString), the more hit strings the 
better. 

Photons originating from farther away from the DOM 
are more likely to have been scattered, and their associ- 
ated distributions of arrival time probabilities are more 
spread out. A larger number of direct hits, that is hits 
that propagate directly to the DOM with little or no 
scattering, has been found to be correlated with better 
track reconstruction. The number of direct hits (NDir) 
is defined as the number of DOMs that have a hit with 
a residual time difference of —15ns < t Tes < 75ns. The 
ratio of direct hits to the total number of detected pulses 
(NDir/NPulses) in an event was also used as a cut pa- 
rameter. 

The length of the event, (LDir), is determined by pro- 
jecting hit DOMs onto the reconstructed track and cal- 
culating the distance between the two endpoints of the 



projection. Larger values indicate a more reliable recon- 
struction of track direction. LDir is calculated using di- 
rect hits only. 

Smoothness is a measure of how well the observed hit 
pattern is explained by the hypothesis of constant light 
emission along the reconstructed muon track. High qual- 
ity tracks have hits equally spaced along the track. This 
parameter, called SmoothAll, is calculated using all hits. 

B. Filtering and Event Selection 

At trigger level, mis-reconstructed atmospheric muon 
events in the zenith region 90° to 180° outnumbered at- 
mospheric neutrinos by a factor of about 10 5 . These mis- 
reconstructed tracks were either individual muon tracks 
or coincident atmospheric muons that mimicked a single 
up-going event. 

Although a variety of filters were deployed at the South 
Pole for the 2008-2009 physics run, events used for this 
analysis were only required to pass the muon filter. The 
muon filter was the primary filter for rejecting down- 
going atmospheric muons, and retaining generic events 
from near or below the horizon. Simple and fast recon- 
structions were performed in real-time at the South Pole. 
These initial reconstructions were less accurate than ones 
performed later, during off-line data processing. How- 
ever, they could be accomplished within the time and 
CPU constraints at the South Pole while keeping up with 
the trigger rate. Zenith angles from LineFit and single- 
iteration SPE likelihood fits, as well as the number of hit 
DOMs (NChannel) and the average number of pulses per 
DOM, were used as selection variables in the muon filter. 
After muon filter event selection was applied, background 
was reduced to a factor of about 10 4 times the neutrino 
event rate. 

Higher level reconstructions included improved likeli- 
hood reconstructions for better angular resolution and 
background rejection, as well as reconstruction of addi- 
tional parameters, such as energy estimation. Fits to 
additional track hypotheses were also performed. Some 
higher level reconstructions incorporated the detailed ice 
model. Prior to higher-level off-line data processing, 
events that were uninteresting or unusable, and that 
clearly were not going to pass final event selection, were 
removed by applying loose cuts based on the results of 
an SPE fit: zenith angle > 80°, RLogL < 12, and PLogL 
< 8. This reduced the amount of background to roughly 
a factor of 10 3 relative to signal. 

Final event selection was accomplished with BDTs that 
used multiple reconstructed and observed parameters as 
input. 

C. Boosted Decision Tree Event Selection 

The Toolkit for Multivariate Data Analysis (TMVA) 
with ROOT [35j was used to implement BDT event clas- 
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FIG. 9: Zenith angles from the SPE fits for the two sub-events created by the geometric split. Neutrino simulation (left) and 
data (right). Box size is proportional to the event density. 



sification. BDTs outperform straight cuts because the 
decision trees are able to split the phase space into a 
large number of hypercubes, each of which is identified 
as either signal-like or background-like 35] . Additionally, 
BDTs often out-perform other multivariate techniques 
because either there are not enough training events avail- 
able for the other classifiers, or the optimal configuration 
(e.g. how many hidden layers for a neural network, which 
variables to use, etc.) is not known and is difficult to de- 
termine [35] ]. Testing with several different multivariate 
algorithms within TMVA indicated that the best results 
for separating signal from background in this case could 
be achieved with BDTs. 

The nodes of a decision tree form what looks like an 
inverted tree. At each node, the algorithm chose the 
particular cut variable and cut value that provided the 
best discrimination between signal and background for 
the events in that node. Events were then split into ad- 
ditional nodes that made up the next layer of the tree, 
and the process repeated until a minimum number of 
events in a node was reached. Variables were used mul- 
tiple times in a tree, with different cut values each time. 
The final nodes were classified as signal or background, 
depending on the classification of the majority of training 
events that ended up in each node. 

Boosting was used to overcome problems associated 
with statistical fluctuations in the simulation used to 
train the BDTs. 200 trees were derived from the same 
training ensemble by re-weighting events. After one tree 
was created, events that were mis-classified in that tree 
had their weights increased, and the next tree was cre- 
ated. This next tree then chose different variables and 
cut values at each node as a result of the altered weights. 
The final classifier used a weighted average of the indi- 



vidual decisions of all 200 trees. 

Two BDTs were used; one having better efficiency 
at lower energies, the other having better efficiency at 
middle and higher energies. Events were accepted if 
their classification score from either BDT exceeded an 
optimized threshold. The function of the BDTs was 
to distinguish between poorly reconstructed background 
events, and signal events that included some that were 
well reconstructed and some that were poorly recon- 
structed. By applying pre-selection cuts prior to training 
the BDTs, some of the poorly reconstructed events were 
removed from the signal event samples, and the overall 
performance of each BDT was improved. For the low 
energy BDT (BDT 1), the pre-selection cut was based 
on LineFit velocity (LineFit velocity > 0.2c). For the 
other BDT (BDT 2), the pre-selection cut was based on 
zenith angles from the split track fits (all four zenith an- 
gles > 80°). The same cuts were applied to the actual 
data as were applied to the simulated background and 
signal event samples used for BDT training and testing. 

Muon neutrino simulation with an E^ 1 spectrum was 
used for signal events in the BDT training. Although the 
true signal spectrum is much steeper than this, testing in- 
dicated this spectrum for training produced a BDT that 
performed better for higher energy events, with no com- 
promise in performance for low energy events. Cosmic 
ray muon simulation from CORSIKA was used for back- 
ground events. Following training, the BDTs were tested 
using independent signal and background event simula- 
tion. Neutrino simulation weighted to an atmospheric 
spectrum, as well as single, double, and triple-coincident 
muon events, weighted to the cosmic ray muon spectrum, 
were used for testing the BDTs. 

Table Q] lists the specific variables used in the BDTs. 
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The NString variable was only used by BDT 1. One ad- 
ditional difference between the two BDTs was the source 
of the Split Track fits. For BDT 1, which was optimized 
for lower energies, the LineFit reconstructions for each 
of the four split tracks (two split geometrically and two 
split in time) were used. For BDT 2, if 16-iteration SPE 
fits were successful for the split tracks, then those results 
were used, otherwise the LineFit results were used. SPE 
fit results were not available for events in which there 
were too few hit DOMs in one or more of the splits to 
perform a likelihood fit. 

TABLE I: Reconstruction variables used in the BDTs. 9z 
refers to the zenith angle. See Sect. IIV"Kl for explanations. 

BDT Variables 

Paraboloid Sigma for the MPE fit 
RLogL from the MPE fit 
PLogL from the MPE fit 
NDir 
LDir 
SmoothAll 
NDir/NPulses 

L °g L Baycs ~ L °g L SPE32 

NString 

|,nMPE flLincFitl 
L °g L Umbr ~ L °g L SPE32 

6z from each of the Split-Tracks 



thetical detector with the same collecting power as Ice- 
Cube, but with 100% efficiency. A c s satisfies the equa- 
tion 

N events = JdtJdSl J dE-<S>{E,e)-A cS {E,6), (5) 

where N events is the number of events passing final se- 
lection cuts and 3>(E,9) is the true flux of atmospheric 
neutrinos with units of GeV" 1 s _1 sr _1 cm~ 2 . In prac- 
tice, the effective area is numerically calculated based 
on the number of neutrino events generated in simula- 
tion, the number passing final event selection cuts, and 
the event weights assigned in simulation to account for 
the probability of reaching and interacting in the detec- 
tor. Figure [TT] shows the effective area as a function of 
energy, for different zenith ranges, at the final cut level. 
Figure IT21 shows the effective areas as a function of energy 
for BDT 1 and BDT 2 separately. 

After eliminating data runs with some strings not op- 
erating, testing in progress, or various faults, there re- 
mained a total of 359 days of livetime from April 2008 to 
May 2009. After final event selection cuts, the number of 
up-going neutrino events was 20,496, with zenith angles 
between 90° and 180°. An apparent excess of horizontal 
events in the data, or deficit in simulation, from 90° to 
97°, will be discussed in more detail in Sect. IV 51 Since 
the origin of this mismatch could not be verified, an addi- 
tional zenith angle cut was applied at 97°. This resulted 
in a data sample of 17,682 atmospheric muon neutrino 
events from 97° to 180°. 



Figure [10] shows the output of each BDT, for the data 
and for simulation weighted to the same livetime (359 
days) . The cut value of 0.73 was chosen to achieve greater 
than 99% purity. Testing the BDTs with simulated signal 
and background data sets indicated that the background 
contamination was less than 0.25%. However, the effec- 
tive livetime of the background simulation available for 
testing was not representative of a year of data. The 
lack of sufficient background simulation near the chosen 
cut values can be seen in Fig. QUI Because we did not 
have a reliable estimate of background contamination, 
comparisons between data and neutrino simulation were 
used to further verify that background rejection was per- 
forming as expected. In particular, the data passing rate 
as a function of BDT cut values was compared to the 
predicted rate from atmospheric muon and neutrino sim- 
ulation. At looser BDT cut values, where sufficient simu- 
lated background events passed the BDT cuts to provide 
a statistically significant estimate, the background from 
simulation underestimated the apparent background in 
the data by about a factor of three. Hence, the amount 
of background contamination in the final data set was 
conservatively estimated to be less than 1%. The ad- 
ditional cut at a zenith angle of 97°, discussed shortly, 
further reduced the potential for background contamina- 
tion. 

The effective area, A e g, is the area occupied by a hypo- 



V. SPECTRUM UNFOLDING 

A. Methodology 

The distribution of the energy-related observable, 
dE/dX, can be expressed as 

b (dE/dX) = A (E v , dE/dX) $ (E v ) , (6) 

where $ is a vector representing the true atmospheric 
neutrino flux as a function of energy, at the point of ori- 
gin in the atmosphere, the vector b is the distribution 
of dE/dX for events in the final sample, and A is the 
response matrix that accounts for the effects of propaga- 
tion through the Earth, interaction in or near the detec- 
tor, detector response, and event selection. An analytical 
solution for A is not known, so it is created from simula- 
tion. 

The desired result from the energy spectrum unfold- 
ing is the true neutrino flux, $. Ideally, this could be 
determined by inverting the response matrix: 

$ (E v ) = A^ 1 (E v ,dE/dX) b (dE/dX) . (7) 

However, direct solution is complicated by the fact that 
events are lost because the detector has limited efficiency 
(many neutrinos either do not interact near the detector 
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FIG. 10: Output of the BDTs for data, as well as neutrino and muon background simulation weighted to the same livetime 
(359 days). The vertical lines mark the chosen cut value of 0.73 for each BDT. 
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FIG. 11: Effective area for up-going muon neutrinos as a func- 
tion of neutrino energy, for various zenith regions. 



or the events do not pass event selection cuts). Addition- 
ally, the detector response is affected by limited energy 
resolution and there is significant smearing of events be- 
tween bins (large off-diagonal elements in the response 
matrix). Moreover, statistical fluctuations in the data 
can lead to unphysical variations in the unfolded spec- 
trum. 

The Sin gula r Value Decomposition (SVD) unfolding 
algorithm [361 ] was used to solve Eqn. [7] and regularize 
the solution. The SVD method involves factoring a non- 
invertible matrix into the product of two orthogonal ma- 
trices and a diagonal matrix, that can then be manipu- 
lated as necessary. This algorithm has been implemented 
in the RooUnfold package [37| for use in the ROOT [H 
data analysis framework. The inputs to the unfolding al- 
gorithm are the response matrix, the predicted histogram 
for the observed distribution, and a histogram for the ex- 
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FIG. 12: Effective areas as a function of energy for each BDT. 
BDT 1 (long-dashed line) performs better than BDT 2 at low 
energies while BDT 2 (short-dashed line) performs better than 
BDT 1 at higher energies. Events are required to pass only 
one of the BDTs and the net effective area is the solid line. 
In contrast to Fig. 111! this plot reflects the zenith-averaged 
effective area for the region 97° to 180°. This corresponds to 
the zenith region used for the analysis. 



pected true flux. 

The expected true flux, $ MC , is a 12 bin histogram 
binned in \og 10 (E^/GeV) from 2 to 5.6, where E u is 
the neutrino energy in GeV. The predicted observables 
histogram, 6 MC , is a 12 bin histogram of the expected 
dE/dX distribution of events passing final cuts, binned 
mlog 10 ((dE/dX)/(GeV/m)) from -2.1 to 1.5. Figurc[]J 
shows the distributions, comparing data to simulation, 
for the observable dE/dX. The response matrix, A, is a 
12 by 12 histogram binned in \og 10 ((dE/dX)/(GeV/m)) 
vs. log 10 (S l/ /GeV), and filled with all events in b MC . 

The response matrix maps the distribution of recon- 
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FIG. 13: Distributions of the dE/dX observable, for data and 
for simulation. The additional cut at zenith angle of 97° has 
been applied. Error bars for data are statistical only. 



structed muon dE/dX values to the distribution of neu- 
trino energies. The correlation between muon energy in 
the detector and the reconstructed dE/dX is not exact, 
and neither is the correlation between muon energy in 
the detector and the incident neutrino energy. The neu- 
trino flux in or near the detector is affected by propa- 
gation through the Earth, during which CC interactions 
attenuate the neutrino flux and NC interactions alter the 
neutrino energy distribution. The muon energy is only a 
fraction of the neutrino energy and only a fraction of the 
muon energy is observed. Below about a TeV, where ion- 
ization rather than stochastic energy losses dominate, the 
energy loss rate is nearly independent of energy. Addi- 
tionally, the stochastic, radiative losses are not uniform 
along the muon's track, as assumed in the reconstruc- 
tion algorithm. If the muon is created in the detector, 
Cherenkov photons generated by the hadronic shower at 
the location of the CC interaction can be detected. If the 
muon is created outside of the detector, it loses some of 
its energy before reaching the detector. 

Figure [14] shows the correlation between neutrino 
energy and reconstructed muon dE/dX. An esti- 
mate of neutrino energy resolution, as a function of 
Log 10 (i? 1 y/GeV), is shown in Fig. [15] To estimate this 
resolution, a Gaussian fit was performed to the dis- 
tribution of Lag 10 ((dE/<2X)/(GeV/m)) in each of 12 
Log 10 (i? ! y/GeV) bins, and the standard deviations from 
these fits are shown in the figure. At higher energies, the 
correlation between neutrino energy and reconstructed 
muon dE/dX is hindered by the fact that the muon 
tracks are not contained within the detector and the 
muon can originate from a significant distance outside of 
the detector. At lower energies, the resolution is aided by 
the fact that the events are more fully contained within 
the detector and the amount of detected light depends 
on the track length within the detector. 

In addition to numerically inverting the response ma- 
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FIG. 14: Correlation between neutrino energy and the recon- 
structed muon dE/dX observable, from simulation weighted 
to the atmospheric neutrino spectrum of [3, 0| • 
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FIG. 15: Estimated neutrino energy resolution, from simula- 
tion. 



trix, the unfolding algorithm applies smoothing (regular- 
ization) to the solution to ensure that statistical fluctu- 
ations in the data do not propagate as unphysical fluc- 
tuations in the unfolded spectrum. The curvature in the 
solution, how sharply it can fluctuate from bin-to-bin, is 
regulated. A regularization parameter enforces a smooth 
cutoff of higher frequency terms in the solution. A lower 
cutoff biases the solution towards the shape of the ex- 
pected spectrum, whereas a higher cutoff allows the so- 
lution to be influenced to a greater extent by fluctuations 
in the data. 

The optimal choice of the regularization parameter de- 
pends on the number of bins and the sample size. Two 
methods for determining the appropriate amount of reg- 
ularization were used, as discussed in Ref. [361 ] - The pri- 
mary method used a result directly from the unfolding 
algorithm, where the coefficients of a particular decompo- 
sition of the rescaled measurement histogram were exam- 
ined. At lower indices these coefficients fall exponentially, 
and the critical term that determines the setting of the 
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FIG. 16: Unfolding of known, toy spectrum. The solid line is 
the assumed spectrum used for regularization, the sum of the 
conventional and prompt atmospheric neutrino flux models 
from Refs. 0,0. The dashed line is the arbitrary, toy spec- 
trum used for generating the toy data. The unfolded result is 
shown without error bars. 



FIG. 17: Results of the atmospheric neutrino energy spectrum 
unfolding. The unfolded spectrum is shown in black, verti- 
cal lines are the estimated uncertainties. The gray line is the 
spectrum that provided the expected shape for the regulariza- 
tion, and includes the conventional atmospheric neutrino flux 
according to Honda et al. 0] and the prompt flux according 
to Enberg et al. [jj. This is the zenith-averaged + flux 
for the region 97° - 180°. 



regularization parameter is at the end of the exponential 
fall, after which the coefficients are not significant. 

As suggested in Ref. [Hj], this result was checked us- 
ing a series of toy simulations that were made system- 
atically and statistically different from the expected true 
distribution. The atmospheric neutrino flux models from 
Refs. [1, Q were used as a baseline. Variations in the 
spectral slope (up to ± 0.1) and normalization about this 
baseline were implemented. For each underlying assumed 
true flux, many randomly fluctuated data sets were gen- 
erated and each simulated data set was unfolded several 
times, using different choices for the regularization term. 
The x 2 °f each unfolded result relative to the true as- 
sumed spectrum was calculated and the distributions ex- 
amined. The regularization term giving the best average 
X 2 was the same as that found by the direct method using 
the decomposition coefficients. 

Fig. [in] shows the performance of the unfolding algo- 
rithm to a toy spectrum. In this example, a toy data set 
was created by arbitrarily modifying the event weights 
in simulation. The spectral slope of the conventional at- 
mospheric neutrino flux (Ref. [3() was made steeper by 
an index correction of -0.05, and the overall normaliza- 
tion was reduced to 80%. Additionally, the prompt flux 
was not included in the toy spectrum, creating a change 
in the shape of the energy spectrum at higher energies, 
where the shape of the actual flux is most uncertain. As 
can be seen from Fig. [T6l some bias is introduced by the 
regularization process at the highest energies where the 
event count is low and the shape of the true spectrum is 
different from the assumed spectrum. 



B. Results and Systematic Uncertainties 

The results of the atmospheric neutrino energy spec- 
trum unfolding can be seen in Fig. 1171 Event selection 
cuts that isolated track-like events, caused by the muons 
created in CC interactions, eliminated localized events 
from the electromagnetic showers induced by v e CC in- 
teractions and the hadronic showers induced by NC in- 
teractions. Additioanlly, production of v T {v T ) by cosmic 
rays is negligible. IceCube is not able to distinguish be- 
tween neutrino and antineutrino events. Hence, the un- 
folded spectrum is the sum of and u^, averaged over 
the zenith region 97° - 180°. 

The major uncertainties in the unfolded spectrum 
are from four categories. These are uncertainties in 
DOM sensitivity and ice properties, zenith-dependent 
data/simulation inconsistencies, statistical uncertainties 
and the impact of the regularization process, and miscel- 
laneous normalization errors such as neutrino cross sec- 
tion and muon energy loss uncertainties. The bin-by-bin 
values for estimates of each of these error sources were 
added in quadrature to obtain the final uncertainty esti- 
mate for each bin of the unfolded flux. 

Systematic uncertainties in ice properties and DOM 
sensitivities lead to systematic errors in the distribution 
of reconstructed dE/dX values, as well as the energy de- 
pendence of the detector's effective area. To estimate 
the impact of these uncertainties, two specialized neu- 
trino simulation datasets were created. In one dataset, 
the number of photons striking each DOM was boosted 
by 10%. In the other dataset, the number of photons was 
reduced by 10%. From this, it was found that a ±10% 
change in the photon flux leads to a ±15% change in 
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event rate and a ±0.09 change in the "apparent" spec- 
tral index of the neutrino flux. These factors were found 
from a three parameter fit that determined the changes 
in normalization, spectral index, and zenith angle tilt, of 
standard atmospheric neutrino simulation, to reproduce 
the best fits to the dE/dX distributions of the simulated 
event samples from these specialized datasets. 

The uncertainty in the DOM sensitivity is taken as 
±8%, based on the measured uncertainty in the PMT 
sensitivity [l6j]. The ±10% change in photon flux in 
the specialized simulation is effectively the same as a 
±10% change in PMT sensitivity, so the normalization 
and spectral index correction factors just mentioned were 
scaled to the ±8% for PMT sensitivity uncertainty. To 
apply these values to ice property uncertainties, the 
change in the average number of photons striking a DOM 
that would result from a change in the ice properties had 
to be estimated. 

First, we assumed a mean propagation length L p = 30 
m with ±10% uncertainty [19j. The propagation length 
for diffusive photon transport from a point source is de- 
fined as 

where L e is the effective scattering length and L a is the 
absorption length. Then, we estimated the fractional 
change in the number of photons at a distance d from 
the union's track as 



where N ~ (l/d)e~ d / Lp is the number of photons 
at distance d for the nominal propagation length, and 
N' ~ (l/d) e - d /( 1± - 1 ) L P i s the number of photons at dis- 
tance d for the perturbed propagation length (nominal 
±10%). The average distance between the track and the 
hit DOMs, per event, was estimated from simulation to 
be about 35 m. The net result of this approximation was 
that the uncertainty in the average photon flux reach- 
ing the DOMs was estimated to be ±12% on average, as 
a result of ice property uncertainties. The normalization 
and spectral index correction factors from the specialized 
simulated data sets were scaled to this ±12% uncertainty 
in the photon flux. 

It should be pointed out that this method of estimat- 
ing the impact of ice property uncertainties is affected 
by two approximations. First, the accuracy of the diffuse 
flux equation is limited at ranges less than several prop- 
agation lengths. Second, changes in ice properties would 
also change the distribution of photon arrival times at 
the hit DOMs, an effect which is not accounted for in 
the specialized simulation data sets. Comparisons were 
made between these simulated data sets and simulation 
generated using Photonics tables derived from a modi- 
fied ice model. In this modified ice model, scattering and 
absorption in the cleaner layers of the ice were arbitrarily 
reduced. This comparison, as well as preliminary results 
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FIG. 18: Possible variability in the true neutrino flux con- 
sistent with DOM sensitivity and ice property uncertain- 
ties. The solid line is the predicted atmospheric neutrino flux 
The dashed lines are the maximum and minimum of 
the possible range of variability consistent with DOM sensi- 
tivity and ice model uncertainties. As mentioned in the text, 
work is on-going to reduce this range of uncertainty. 



from ongoing work to improve the simulation of photon 
propagation in the ice and derive a more accurate esti- 
mate of uncertainties related to ice properties, indicated 
that the method used here likely over-estimates the im- 
pact of ice property uncertainties on the normalization 
and apparent spectral index, particularly in the higher 
energy bins. 

Adding the uncertainties in detector response due to 
ice properties and DOM sensitivity in quadrature leads 
to an estimated ±22% uncertainty in the normalization, 
correlated with an uncertainty of ±0.13 in the apparent 
spectral index. These detector uncertainties lead to un- 
certainties in the apparent neutrino flux. For a given 
detector response, i.e. a measurement of the dE/dX dis- 
tribution, the true normalization and spectral index of 
the neutrino flux cannot be constrained better than al- 
lowed by these uncertainties. Figure HH] shows the re- 
sulting range of uncertainty in the measurement of the 
atmospheric neutrino energy spectrum. 

As mentioned in Sect. IIV CI there was a statistically 
significant excess of events in data, or a deficit in simula- 
tion, between 90° — 97°, i.e near the horizon. Figure IT51 
shows the cos(zenith) distributions for data and for simu- 
lation, with simulation normalized to the data. A similar 
excess was also observed in the AMANDA detector [39| . 
A number of checks and tests were performed, includ- 
ing evaluation of track quality parameters, the depth- 
dependence of the excess, the strength of the BDT scores, 
and visual examination of a subset of events in a software- 
based event viewer. The horizontal excess in data does 
not decrease with depth, nor with tightened BDT cuts. If 
the BDT cut is loosened, mis-reconstructed muons show 
up predominantly near the top of the detector, as ex- 
pected. These checks are consistent with the possibil- 
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FIG. 19: Cosine(#z) distributions for data and for simulation, 
using zenith angle from the MPE fit. Simulation has been 
normalized to the data. Error bars for data are statistical 
only. 



ity that the excess events are due to muons from atmo- 
spheric neutrino interactions. However, it is also possible 
that they are due to an excess of mis-reconstructed at- 
mospheric muons. 

It is likely that the lower event rate in simulation, close 
to the horizon, is due to uncertainties in the simulation 
of Cherenkov photon propagation in the ice or of inaccu- 
racies in the simulation of cosmic ray events, such as in- 
sufficient livetime, limitations with the cosmic ray model 
or its implementation in CORSIKA, or uncertainties in 
muon propagation and energy loss. Hence, it cannot 
be excluded that the horizontal excess is due to resid- 
ual and un-simulated atmospheric muons and coincident 
events. It could also be related to uncertainties in the 
atmospheric neutrino flux due to atmospheric variability, 
discussed shortly. Since we were not able to verify the 
precise origin of the mismatch near the horizon, events 
in the zenith region 90° to 97° were not used. 

To estimate the impact of any remaining zenith- 
dependent systematic uncertainties in the zenith range 
97° — 180°, separate unfoldings were performed for the 
zenith range 97° - 124° and the zenith range 124° - 180°. 
The results of this test are shown in Fig. [20l together 
with the predicted zenith-averaged flux corresponding to 
each angular range. The differences between result and 
prediction are not consistent between the two regions. 
For the more vertical events (gray in Fig. [20]) . the flux is 
lower than predicted for middle and higher energies. For 
the more horizontal events (black in Fig. |2"U|) . the flux is 
slightly lower than predicted at low and at high energies. 
The relative differences between result and prediction for 
the two zenith regions was taken as an estimate of the 
impact of anisotropic uncertainties. 

Seasonal and regional variations in the atmospheric 
temperature profile are expected to lead to variations in 
the atmospheric neutrino flux [40j], and could be caus- 




Log io (E v /GeV) 



FIG. 20: Comparison of unfolded energy spectra for differ- 
ent zenith ranges. Separate unfoldings were performed for 
the zenith range 97° — 124° (black) and the zenith range 
124° — 180° (gray). The unfolded results for each region (hor- 
izontal lines) and the predicted spectrum corresponding to 
each region (curves) are shown. Uncertainties for these re- 
sults are not shown. 

ing the direction-dependent differences between data and 
simulation. Colder temperatures correspond to a greater 
air density and a shallower atmosphere. Greater at- 
mospheric density leads to more collisions of pions and 
kaons prior to their decay. Hence, the production of high 
energy neutrinos is reduced. The converse occurs for 
warmer temperatures. The kinematics of collision and 
decay, and slant angle through the atmosphere, conspire 
to lead to variations in the energy and zenith angle depen- 
dence of atmospheric neutrino production for different 
atmospheric conditions. The normalization uncertainty 
on the Honda et al. conventional atmospheric neutrino 
flux model includes an estimated 3% uncertainty due to 
uncertainties in the atmospheric density profile [!, l4lj . 
However, the flux calculation uses a climatological av- 
erage atmosphere (the US-standard '76). The estimate 
of the error in the flux calculation is based on the error 
in the climatological average atmospheric density profile. 
It does not account for changes in normalization, or in 
the energy and zenith distribution of atmospheric neutri- 
nos, that results from regional and seasonal atmospheric 
variability. 

The impact of statistical uncertainties in the data, as 
well as bias due to the regularization process and the pos- 
sibility that the assumed spectrum used to compute the 
amount of regularization may be different from the true 
spectrum, were estimated using toy simulations. First, 
a six-parameter forward folding fit to the data was per- 
formed. In the forward folding fit, the general form of 
the flux was assumed to be consistent with the shape of 
the theoretical predictions 0, , but corrections to the 
normalization, spectral index, and zenith angle tilt of the 
conventional and prompt atmospheric neutrino flux mod- 
els were propagated through simulation. The fit variables 
that produced the best fit between the simulated detector 
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FIG. 21: Statistical and regularization-induced uncertainties 
in the unfolded result. The errors in each bin are given as the 
percent of the true flux. 



response and the data were used to re-weight simulated 
events in the toy experiments to mimic the data. The 
results of the forward folding fit indicated a possible sys- 
tematic suppression of the neutrino event rate at higher 
energies, and this suppression was included in the toy 
simulations. 

One thousand trials were performed, with events in 
each bin of the toy dE/dX distributions fluctuated ac- 
cording to a Poisson distribution. Statistical uncertain- 
ties in the neutrino simulation were also included. The 
difference between the unfolded energy spectrum and the 
known, "true" spectrum that the toy experiments were 
based on were computed for each trial. The 68 th per- 
centile of the errors in each bin from the 1000 trials were 
assigned as the uncertainty. The result of this analysis 
of statistical and regularization uncertainties is shown in 
Fig. 1211 where the errors in each bin are given as percent 
of the true flux. A potential systematic bias between 
the shape of the true flux and the shape of the assumed 
flux used to train the unfolding algorithm accounts for 
roughly half of the uncertainty indicated in Fig. [5T] for 
the two highest energy bins. 

A 3% uncertainty in the charged current, deep-inelastic 
neutrino-nucleon scattering cross section is estimated to 
lead to a 3% uncertainty in atmospheric neutrino event 
rates, and uncertainties in muon energy loss are esti- 
mated to lead to a 1% uncertainty (42[. Reconstruction 
and cut biases are estimated to introduce a 2% uncer- 
tainty in event rate. Adding these, and the 1% back- 
ground contamination, in quadrature gives a 4% uncer- 
tainty in the event rate, assumed to be independent of 
energy. 

A summary of uncertainties in the unfolded result can 
be seen in Fig. [251 as wen as Table [TTJ At the lower 
end of the unfolded energy range, uncertainties are domi- 
nated by zenith-dependent inconsistencies. At the middle 
of the range, uncertainties are dominated by the DOM 
sensitivity and ice property uncertainties, as well as the 



FIG. 22: Sources of uncertainty in the unfolded energy 
spectrum. The solid lines are the systematic uncertainties 
due to DOM sensitivity and ice property uncertainties, the 
short-dashed lines are the uncertainties implied by zenith- 
dependent inconsistencies in data/simulation comparisons, 
and the long-dashed lines are the statistical and regulariza- 
tion uncertainties from toy MC studies. Not shown is the 
uniform 4% uncertainty due to miscellaneous normalization 
errors assumed to be independent of energy. 

zenith-dependent uncertainties. Uncertainties in DOM 
sensitivity and ice properties dominate at higher ener- 
gies. 



TABLE II: Zenith-averaged, unfolded atmospheric muon neu- 
trino energy spectrum. 



log 10 (Ev/GeV) 


dN/dE v ■ El 
(GeVs^sr^cm" 2 ) 


% Uncertainty 


2.0 - 


2.3 


3.6 x 10" 4 


+29, -28 


2.3 - 


2.6 


1.6 x 10" 4 


+21, -22 


2.6 - 


2.9 


7.0 x 10" 5 


+31, -32 


2.9 - 


3.2 


2.8 x 10~ 5 


+50, -50 


3.2 - 


3.5 


1.1 x 10" 5 


+65, -62 


3.5 - 


3.8 


4.0 x 10~ 6 


+71, -63 


3.8 - 


4.1 


1.4 x 10" 6 


+74, -58 


4.1 - 


4.4 


4.7 x 10~ 7 


+82, -53 


4.4 - 


4.7 


1.6 x 10" 7 


+95, -53 


4.7 - 


5.0 


5.4 x 10 -8 


+113, -57 


5.0 - 


5.3 


2.0 x 10~ 8 


+135, -64 


5.3 - 


5.6 


7.9 x 10~ 9 


+158, -72 



VI. CONCLUSIONS AND OUTLOOK 

A zenith averaged unfolding of the atmospheric muon 
neutrino flux {v^ plus v^), from 100 GeV to 400 TeV, 
was performed. This is the first atmospheric neutrino 
measurement to such high energies, and the spectrum 
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FIG. 23: Comparison with previous measurements of the atmospheric neutrino energy spectrum; the Frejus result [43fl. upper 
and lower bands from SuperK an AMANDA forward folding analysis [45l ]. and an AMANDA unfolding analysis [39t] . All 
measurements include the sum of neutrinos and antineutrinos. The AMANDA unfolding analysis was a measurement of the 
zenith-averaged flux from 100° to 180°. The present analysis (IC40 unfolding), which extends the measurement up to 400 TeV, 
is a measurement of the zenith-averaged flux from 97° to 180°. Vertical error bars include systematic as well as statistical 
uncertainty. 
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Log io (E/GeV) 

FIG. 24: Comparison of various prompt flux models to the 
unfolded result. The models shown are the sum of the Honda 
flux Q], plus one of Sarcevic 0, Naumov [llj]. or Martin [Iol |. 
Vertical error bars include systematic as well as statistical 
uncertainty. 



is consistent with predictions for the atmospheric muon 
neutrino flux. However systematic uncertainties will need 
to be reduced before specific flux models [D, [U, l9Hll1| can 
be constrained. In particular, we are as yet unable to 
confirm the contribution of a prompt flux. Figure [53] 
compares the results of this analysis (IC40 unfolding) to 
previous measurements of the atmospheric neutrino en- 
ergy spectrum. As discussed in Sect. IV Bl the estimate 
of uncertainties in the IceCube result are dominated by 
DOM sensitivity and ice property uncertainties, as well 
as the zenith-dependent mismath between data and sim- 
ulation. These uncertainty estimates are expected to be 
reduced as our simulation is improved. A comparison 
between the unfolded spectrum and various prompt flux 
models @-[H| is shown in Fig. [24l 

Several improvements are anticipated in atmospheric 
neutrino measurements with IceCube. Correlations be- 
tween variations in atmospheric temperature profiles, 
and the energy and zenith angle dependence of the atmo- 
spheric neutrino flux, are being investigated using in-situ 
atmospheric temperature measurements. Pulsed LED 
sources installed on each DOM are being used to ex- 
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tend the ice description to the deepest ice in the detector 
with in-situ measurements like those done in AMANDA 
for the ice down to 2100 m. Studies with cosmic ray 
muons are being used to reduce the uncertainty in DOM 
sensitivity. Work is also on-going to identify and cor- 
rect potential problems in simulation that could be con- 
tributing to data/simulation mismatch. Perhaps most 
significantly, this includes improving the simulation of 
light propagation within the detector, which is antici- 
pated to improve the data to simulation agreement for 
several measured and reconstructed variables. These im- 
provements will be discussed in a future paper. Once 
simulation of light propagation in the ice has been im- 
proved, it should be possible to use a more sophisticated 
and realistic method for estimating the impact of ice 
model uncertainties. As data collection continues, and 
improvements to simulation are implemented, it will be 
possible to extend the measurement of the atmospheric 
neutrino energy spectrum with IceCube to PeV energies, 
as well as significantly reduce the uncertainties. 
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